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1 Introduction 



Synchronization has recently been attracting many researchers' attention as a 
simple yet interesting example of collective behavior on complex networks [1- 
22]. Of particular interest is how the networks' ability to synchronize depends 
on various structural parameters of the networks, such as average node-to- 
node distance [23,24], clustering coefficient [25], degree distribution [24,26], 
and weight distribution [27-29]. Revealing the precise mechanism for syn- 
chronization in relation to the network structure is an important step toward 
understanding more complex behavior and unveiling the evolutionary origin 
of real- world networks. 

In some naturally evolved networks, such as neuronal and biochemical net- 
works, there is evidence that synchronized or more general coordinated be- 
havior may be playing signicant roles in the system's functionality [30-32]. 
It appears natural to expect that the ability of such networks to synchronize 
their activity has been optimized to some extent during their evolution. In re- 
ality, however, the problem is likely to be more complicated, since the fitness 
of the networks could depend on multiple factors, such as stability, robustness, 
and adaptability. Because of the complexity of the problem, researchers have 
instead been focusing on simplified yet tractable optimization problems as a 
first step toward solving the real problem. In particular, there have been some 
efforts to solve the problem of maximizing the synchronizability of oscillator 
networks [9,21,22,27-29,33]. The scope of such studies has so far been mainly 
limited either to numerical investigations or to single-parameter families of 
possible networks. 

In this paper, we go beyond these restrictions and present rigorous solutions 
to the problem of maximizing the network synchronizability, measured by the 
range of coupling parameter for which the system achieves stable synchro- 
nization. We also consider a more general concept, the cost required for stable 
synchronization, and treat the problem of minimizing it. Remarkably, we prove 
that the solution sets of the two optimization problems coincide and are com- 
pletely characterized by a simple condition on the Laplacian eigenvalues of 
the network. This spectral characterization, however, does not provide much 
intuition about the structure of the optimal networks. To gain more intuition, 
we explicitly construct a large subclass of optimal networks characterized by a 
hierarchical structure, in which information can fiow only from top to bottom 
of the hierarchy, making the network optimal for synchronization. 

The theorems in this paper also provide mathematical foundation for the so- 
lutions of a related problem that we presented in our recent publication [34]. 
That problem was originally motivated by the discovery that random scale- 
free and other degree-heterogenous networks are generally difficult to syn- 
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chronize [24]. This discovery led to efforts to enliance the synchronizabihty of 
complex networks by introducing directionality and weight to each link [27- 
29]. Underlying such efforts is a problem of optimization with topological con- 
straints [34] : given a fixed topology of allowed interactions, find assignments of 
weights and directions that would maximize the range of the coupling param- 
eter for stable synchronization. Here we prove that a solution can be system- 
atically and explicitly constructed using oriented spanning trees embedded 
within any given connected topology of allowed interactions. The resulting 
networks are guaranteed to be optimal with respect to both synchronizabihty 
and synchronization cost. 

The problem of optimization under topological constraints is potentially rel- 
evant for many real-world networks. In metabolic networks, for example, the 
weights and directions of feasible links (metabolic fluxes) are adjusted to op- 
timize fitness, which is likely to account for robustness of synchronized be- 
havior against environmental changes [35]. Similar adjustment of weights and 
directions may enhance neuronal synchronization within a given topology of 
synaptic connections in the brain. In designing the interaction scheme for a 
computer networks, choosing proper weights and directions may optimize the 
performance of computational tasks based on synchronization of processes [36]. 
The adjustment of fiows in power grids and communication patterns in social 
organizations are additional examples where directional and weighted patterns 
may be favored because they can facilitate the synchronized or coordinated 
behavior on which the functioning of these networks is based. 

Our results are based on an extension [34] of a well-known framework for 
studying network synchronization [37]. The power of this framework is that it 
can separate the effect of the network structure from that of the dynamics of 
individual nodes. However, it implicitly assumes that the Laplacian matrix of 
the network is diagonalizable, i.e., the dynamics of the network must be de- 
composable into independent eigenmodes. In most of the previous works, this 
assumption was automatically satisfied, since the main focus was on symmet- 
rically coupled networks, which are guaranteed to be diagonalizable. However, 
the same does not hold true in general when the network is directed. In fact, we 
prove the interesting result that most optimal networks are non-diagonalizable, 
and thus violate an assumption of the original framework. We show that the 
stability condition for synchronization is formally the same for all networks, 
but the speed at which the system converges toward the synchronized state 
can be significantly slower when the network is non-diagonalizable. 

The technique underlying our extended framework can be regarded as an 
example of a methodology for studying complex systems that relies neither on 
the eigenmode decomposition nor on any kind of superposition principle, and 
is expected to meet applications in various other network phenomena. 
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The paper is organized as follows. In the next section, we present the extended 
master stability framework. We state the optimization problems in Section 3 
and characterize their solutions in Section 4, 5, and 6. In Section 7, we consider 
the optimization problems with topological constraints. Finally, in Section 8, 
we make concluding remarks on implications of our results and on future 
directions. 



2 Extension of Meister Stability Analysis 



Consider n identical oscillators whose individual dynamics without coupling is 
governed by x = F(x), x e IRJ^. Now consider the network of these oscillators 
interacting through a diffusive-type coupling, i.e., oscillator i receives input 
from oscillator j that is proportional to Ay[H(xj) — H(xj)], where Aij is a 
nonnegative constant representing the relative strength of the coupling, and 
H : IRJ^ — > IRJ^ is a general output function. The interaction is indeed the 
usual diffusive couphng when H(x) = x. The set of equations governing the 
dynamics of the system is then 

n 

Xi = F(xi) + (7^Ai,[H(x,-)-H(xi)], i = (1) 

where a is the parameter controlling the overall coupling strength. Note that 
the system can be regarded as a sum of two distinct components: the network 
structure represented by the adjacency matrix A — (Aij) and the dynamical 
component represented by the functions F and H. The general method of anal- 
ysis introduced by Pecora and Carroll [37] to study the stability of complete 
synchronization in Eq. (1) is based on the diagonalization of its variational 
equation. In the following, we extend their analysis to include cases where the 
diagonalization is not necessarily possible. 

We first note that having diffusive-type coupling guarantees the existence of 
a completely synchronous state, though it may be unstable. In fact, given 
any solution x = s{t) of the individual dynamics x = F(x), the completely 
synchronized solution, defined by Xj = s{t), i — is automatically 

a solution of the entire system (1). For notational convenience, we rewrite 
Eq. (1) as 

n 

x, = F(x,)-a^L,,H(x,), (2) 
where L — (Lij) is called the Laplacian matrix of the directed weighted net- 
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work, defined by 



Lij = -Aij ifi^j, 

La — — ^ Lij. (^) 

Like the adjacency matrix A, the Laplacian matrix also contains all informa- 
tion about the network structure and can be regarded as a network analog of 
the Laplacian operator for diffusive processes on continuous space. Note that 
L is not necessarily symmetric because in our general formulation the network 
is not constrained to be undirected. 

The stability condition can be studied by considering the variational equation 
for the synchronous solution Xj = s{t) of Eq. (2): 

n 

ii = ^(s)^, - a ^ LijDH{s)Cj, (4) 
which can also be written in the matrix form as 



(5) 



where ^ = i^i,---,^n) m x n perturbation matrix, is the vector 

of perturbations to the ith oscillator, and L^ denotes the transpose of L. In 
the Pecora-Carroll analysis, the assumption that the Laplacian matrix L is 
diagonalizable was implicitly used to diagonalize the variational equation (5). 
Here we do not assume the diagonalizability of L. Instead, we utilize the 
Jordan canonical transformation of L. For any nxn matrix L, there exists an 
invertible matrix P of generalized eigenvectors of L, which transforms L into 
Jordan canonical form as 



P-^LP = J 







Bi 



\ 



B, 



(6) 



where S, 's are blocks of the form 



Bi 



1 A 



(7) 



1 A 
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and A is one of the (possibly complex) eigenvalues of L. We note that the 
Jordan canonical transformation has been used to study the stability of syn- 
chronization in specific classes of networked systems [38-40]. By applying the 
change of variable r] = i{P~^Y Eq. (5), we get [34] 

7) = DY{s)r] - (7DH(s)77 J^. (8) 

Each column of r/, being a linear combination of all ^^'s, represents in general a 
mode of perturbation to the entire oscillator network, and not to any particular 
oscillator. Thus, the synchronous solution is stable if and only if all of these 
columns converge to zero under Eq. (8). 

Before getting into the general treatment, let us first consider the case where 
L is diagonalizable. In this case, the matrix J is a diagonal matrix having the 
eigenvalues Ai, . . . , A„ of L along the diagonal. Thus, the equation for each 
column of r] becomes independent of the others and takes the form 

y = [DF(s) - aL'H(s)]y, (9) 

where a — crAj when y represents the ith column of rj. Regarding a as a tun- 
able complex parameter, Eq. (9) is called a master stability equation and its 
stability profile as a function of a determines the linear stability of the syn- 
chronous solution in systems with various network structures represented by 
the Laplacian eigenvalues. The largest Lyapunov exponent A(q;) for the solu- 
tion y = is usually used to test the stability and is called a master stability 
function [37]. The eigenvalues Ai, . . . , A„ of L can always be arranged so that 
= Ai < Re A2 < ... < Re An, since Lij = implies that we always have 
eigenvalue corresponding to the eigenvector (1, . . . , 1)^, and all eigenvalues 
are guaranteed to have nonnegative real parts by the Gerschgorin Circle The- 
orem. Thus, the condition for the synchronized solution to be linearly stable 
is 

A(aA,) < 0, i = 2,3,...,n. (10) 

See Fig. 1 for a visual demonstration of this stability condition. Note that Ai = 
is excluded from the condition, because A((tAi) = A(0) actually determines 
the linear stability of the individual solution s(t) against perturbation [try 
setting q; = in Eq. (9)]. It would be positive if s(i) is chaotic, but it does not 
affect the stability of synchronization. In other words, A(0) corresponds to the 
stability in the direction parallel to the synchronization manifold defined by 
{xi = • • • = x„}, while A(crA2), . . . , A(crA„) correspond to the stability in the 
directions transversal to the manifold. 

Let us now consider the more general case where L is not necessarily diag- 
onalizable. Each block of the Jordan canonical form corresponds to a subset 
of these columns in r), which obeys a subset of equations in (8). For example, 
if block Bi is k X k, and if the corresponding columns of 77 are denoted by 



6 








Large o 







Re a 



Re a 



Re a 



Fig. 1. Schematic illustration of the stability condition (10) for synchronized state. 

The shaded area is the region in the complex plane in which the master stability 
function A(a) is negative, and the dots represent aXi for i = 2, . . . , n. The condi- 
tion (10) corresponds to having all the dots in the shaded region. 



T]i,r]2, . . . ,rij^, then the equations take the form 

77i = [DF{s) - aDIl{s)]rj, 

= [DF{s) - aDH{s)]ri2 - aDH{s)rj, 



(11) 
(12) 



f], = [DF{s) - aL>H(s)]77, - aDIl{s)rj,_,, 



(13) 



where a = aX. Here 'ni-,'n2-> ■ ■ ■>''lk represent the modes of perturbation in 
the generalized eigenspace associated with eigenvalue A. Equation (11) has 
exactly the same form as the master stability equation (9), so r)i converges 
exponentially to zero as i — > oo if and only if A((tA) < 0. The condition for 
Eq. (12) to be stable is apparently more involved but can be formulated as 
follows. Assuming that A{aX) < and that the norm of DH(s) is bounded, we 
have that the second term in Eq. (12) is exponentially small as well. Then, the 
same condition A{aX) < guarantees the stabihzing effect of both the first and 
second terms, resulting in exponential convergence of 772 to zero as t ^ 00. The 
same argument applied repeatedly shows that 773, . . . , 77^ must also converge 
to zero if A{aX) < 0. This shows that A{aX) < is the condition for the 
linear stabihty of the equations corresponding to each full block Bi. Thus, 
when all the Jordan blocks are taken into account, we see that the stability 
condition for the synchronous solution in the general non-diagonalizable case 
is also given by (10). 

Although the stability condition is the same for both the diagonalizable and 
non-diagonalizable cases, it is worthwhile noting that there is a crucial dif- 
ference. If L is diagonalizable, then each mode of perturbation is decoupled 
from others, so the exponential convergence of each column of rj occurs simul- 
taneously and independently of other columns. On the other hand, if L is not 
diagonalizable, some modes of perturbation may suffer from a long transient 
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because they may be coupled to other modes associated with the same eigen- 
value. In Eqs. (11)-(13), may have to wait for rj^ to get small enough before 
it can start converging, 773 may have to wait for 772, and so on, so 77^ may have 
to wait for a long time before it starts to converge. The larger the size k of the 
Jordan block, the longer we expect the transient to be. As a simple example 
to illustrate this effect, consider a network of linearly coupled phase oscillators 
described by 



UJ 



i=i 



uj-a^LijOj, 



(14) 
(15) 



where §^ denotes the unit circle. In this case, the corresponding variational 
equation around the synchronized solution 6i = ujt, i = 1, . . . ,n, is a simple 
linear system, and so are the corresponding Eqs. (11)-(13). Thus, they can be 
explicitly solved to give 



V3 



Vk 



cie 



-at 



-at 



-ciat + C2)e 

C2(Tt + C3 



-at 



-1) 



{k-l)\ 



+ 



-at 



(16) 
(17) 

(18) 

(19) 
(20) 



where (ci, . . . , c^) is the initial value of the perturbation vector (t^i, . . . , 77^). 
We can indeed see that the polynomial factors lead to slower convergence for 
larger k in this example. 

In fact, even when L is diagonalizable, we expect to see longer transient as 
it becomes closer to being non-diagonalizable, i.e., some of its eigenvectors 
become closer to being parallel. To see this, imagine a small sphere centered at 
the origin of the space of perturbation matrices ^. If a pair of eigenvectors are 
almost parallel, then the matrix P of eigenvectors is close to being singular. 
Hence, the sphere is stretched quite a bit along some direction under the 
transformation 7] = ^{P"^)^. This implies that small perturbations to the 
synchronized state in the original coordinates can lead to large perturbations in 
the eigenvector coordinates. This in turn leads to relatively long transient, even 
though the type of convergence remains purely exponential. This mechanism 
is at work to some extent for any network with non-orthogonal eigenvectors, 
but the effect is more prominent if the eigenvectors are closer to being parallel. 
In the limit of parallel eigenvectors, L becomes non-diagonalizable, and the 
convergence becomes qualitatively different, as we saw in the phase oscillator 
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example above. Thus, wc expect to observe relatively long transient not only 
in a few special networks, but also in many others close to them. 

The fact that the stabihty condition is the same regardless of the diagonahz- 
ability of the Laplacian matrix is analogous to the fact that the linear stability 
condition of a fixed point in a dynamical system is the same regardless of the 
diagonalizability of the Jacobian. In both situations, lack of diagonalizability 
leads to long transient. 



3 Optimization Problems 

Here we are interested in two different optimization problems: 

• Which network structure maximizes the synchronizability of the system? 

• Which network structure allows the system to synchronize stably with the 
minimum possible cost? 

To address these optimization problems, we need to precisely define quan- 
tities to optimize: the synchronizability and the synchronization cost of an 
oscillator network. To set the stage for doing this, we first let TZ denote the 
stability region^ defined as the subset of the complex plane in which the master 
stability function A(a) is negative, i.e., 7^ = {« G C | K{a) < 0}. Using this 
notation, the stability condition (10) can now be written as aA2, . . . , crA^ G 7?.. 
This shows clearly that there are only two distinct factors that determine the 
stability of the synchronized solution: 

(1) Network structure, encoded in the adjacency matrix A and affecting the 
stability only through the Laplacian eigenvalues A2, . . . , A„; 

(2) Dynamical component, consisting of the individual dynamics given by 
F and s, together with the output signal function H, and affecting the 
stability only through the properties of the stability region TZ. We denote 
this component by (F,H,s). 

By considering the complex conjugate of Eq. (9), we can see that TZ is always 
symmetric about the real axis. In most of the previously studied cases, it has 
been found [37,41] that 

(A-I) 7^ is a convex subset of C. 

In addition, for a large class of systems in which the dynamics of each oscillator 
is chaotic, it has also been found [18,28,37,41] that 

(A-II) < CKi < 0:2 < 00, 
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where we define 



Q!i = inf{Rea | a G 71} , a.^ = sup{ReQ; | a G TZ}. (21) 

The assumptions (A-I) and (A-II) on the stability region imply that the set 
of overall coupling strength a for which complete synchronization is stable, 

hync = {o- 1 A((7Ai) < 0, i = 2, . . . , n} 
= {a I aX2, aXn G 7^}, 

is either an empty set or a finite interval with endpoints at CTmin and o"max) where 
< cxmin < cTmax < cxD. lu the case of a finite interval, this can be physically 
interpreted as follows. The completely synchronized state of the network is 
unstable for small enough values of a, but as a is increased, it becomes stable 
at a lower threshold a^^i^ and then becomes unstable again above an upper 
threshold (Tmax- Thus, the relative width of this interval defined by 

(23) 

'-'^min 

provides a natural and convenient measure of how easy it is for the network 
to synchronize, i.e., the synchronizability of a network of coupled chaotic os- 
cillators: larger values of S correspond to more synchronizable networks. In 
this paper, we consider only those systems with TZ satisfying properties (A-I) 
and (A-II) to ensure that S is well-defined. 

The synchronization cost C of a network is defined [27,28] as the sum of the 
total input strength of all nodes at the lower synchronization threshold amin' 

n 

C = a„,in Yl ^ij- (24) 

We define 5 = and C = oo when Isync is empty because this simply means 
that stable synchronization is impossible. Note that while S is guaranteed to 
give meaningful values only when (A-I) and (A-II) are satisfied, C is mean- 
ingful without any assumptions on the stability region. This means that the 
notion of synchronization cost applies to a wider range of systems, includ- 
ing those with no short-wavelength bifurcation [19,42] or with intermediate- 
wavelength bifurcation. See Fig. 2 for examples of the stability region for such 
systems. The results presented below for the synchronization cost would hold 
even when these assumptions are weakened, but here we assume (A-I) and 
(A-II) for simplicity. 

With S and C defined, the two optimization problems under consideration 
can now be formulated precisely: 
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(a) 



(b) 



Ima 



Ima 








Re a 



Fig. 2. Examples of stability region TZ for systems with (a) no short-wavelength 
bifurcation or with (b) intermediate-wavelength bifurcation. 

• Which network structure A has the property that it maximizes the synchro- 
nizability S for any given dynamical component (F, H, s) with TZ satisfying 

(A-I) and (A-II)? 

• Which network structure A has the property that it minimizes the syn- 
chronization cost C for any given dynamical component (F, H, s) with TZ 
satisfying (A-I) and (A-II)? 

By requiring the optimality for the entire class of dynamical components, 
we are defining the optimality of a network structure, independently of the 
dynamical component of any specific system. Note that both S and C are 
invariant under re-scaling of A, and thus only the relative distribution of the 
individual coupling strength Aij is important for the optimization problems. 



4 Optimal Network Structures 

In this section, we give a complete characterization of the class of solutions to 
the two optimization problems introduced in the previous section. We start 
by studying the properties of the quantities to be optimized. Even though the 
synchronizability S and the synchronization cost C generally depend on both 
the dynamical component encoded in TZ and the network structure encoded in 
A2, . . . , A„, it turns out that they are both bounded by a quantity that depends 
only on TZ (and also on n in the case of C) . 

Theorem 1 For any stability region TZ satisfying (A-II), we have 



PROOF. Let the Laplacian eigenvalues be ordered so that = Ai < Re A2 < 
• • • < Re A„. If Re A2 = 0, then /gync is empty, and the inequalities are clearly 



S<—, C>a^{n-1). 



(25) 
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satisfied since S = and C = oo. Thus, we are left with proving the inequah- 
ties in the case < Re A2 < • • • < Re A„. 

Suppose a e /sync- Then, since aXi & TZ ior i — 2, . . . , n, we have 

ai < Re((7Ai) < q;2, (26) 

and hence 

<a<-p- (27) 



Re Xi Re Aj 

for i — 2, . . . ,n. In particular, this implies that 



<<7<-^. (28) 



ReA2 Re A 

This holds for any a e /sync, so by the definition of CTmin and Umax, we have 



< CTmin and (Tmax < TT^" (^9) 



Re A2 Re A„ 



Therefore, 



and 



< ^ . ^ < ^, (30) 
Q!i Re A„ cti 



C — Cjjjjjj — CTmin • tr L 



= CTmin XI - p X ^® ''^^ 

i=2 ^6 ^2 j=2 

> -^•(n-l)-ReA2 = Q;i(n-l). □ (31) 
Re Ao 



In the special cases where all the Laplacian eigenvalues are real, 5" and C 
can be computed explicitly. In particular, S is proportional to an eigenvalue 
ratio [23]. Such situation occurs for example for the undirected networks, for 
which the Laplacian matrix is symmetric. This result is expressed in the fol- 
lowing theorem. 

Theorem 2 Suppose that the Laplaeian eigenvalues are all real and ordered 
as = Ai < A2 < • • • < A„. For any TZ satisfying (A-I) and (A-II), we have 

S^^-^. C^^±K. (32) 

Oil ■^n ^^2 i=2 



PROOF. Suppose TZ satisfies (A-I) and (A-II). Let T^reai be the intersection 
of TZ and the real axis. By combining (A-I) with the fact that TZ is symmetric 
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about the real axis, we see that TZreai is an interval with endpoints at ai and 
a2- We can write 

n 

/sync = n 4?nc, where = {a I aX, G TZ}. (33) 

Since Aj is real, we have crAj e 7?. if and only if aAj is in the interval 7?.reab ^-iid 
hence /g^y^^^, is an interval whose endpoints are at cci/Aj and a2/Xi- Taking into 
account all i = 2, . . . , n, this means that /gync is an interval with endpoints at 
c^min = ai/X2 and a^ax = Qi2/A„. Thus, 

5 = ^ = ^ . ^. (34) 

Cmin Oil A„ 



We also have 



C = a„,inE A, = ^-trL=^X^A,. □ (35) 



A surprising consequence of this theorem is that a simple condition on the 
eigenvalues suffices to guarantee both the maximum synchronizability and the 
minimum synchronization cost for any dynamical component (F, H, s) with 
TZ satisfying (A-I) and (A-II): 

Corollciry 3 Suppose that the Laplacian eigenvalues of a network satisfy, 

= Ai < A2 = • • • = A„. (36) 
Then, S and C achieve their maximum and minimum values, respectively, i.e., 

S^—, C = ai{n-1), (37) 

Oil 

for any (F, H, s) with TZ satisfying (A-1) and (A-11). 



PROOF. Condition (36) implies that the eigenvalues are all real, so Theo- 
rem 2 applies, and we have 

S = — ■ = —, (38) 

«! Xn 



and 



(:7=-i^A, = ai(n-l). □ (39) 

^^2 i=2 



Corollary 3 shows that any network satisfying (36) is a simultaneous solution 
of the two optimization problems under consideration. However, even more 
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surprising is the fact that those that satisfy (36) arc actually the only solutions. 
In other words, the two classes of optimal networks — those with the maximum 
synchronizability and those with the minimum synchronization cost — actually 
coincide, and both can be completely characterized by condition (36). To state 
this in a precise but convenient form, let us define the following terminology. 
We say that a network given by A has the maximum synchronizability ii S — 
a2/<yi for any (F,H,s) with TZ satisfying (A-1) and (A-11). Similarly, we say 
that a network given by A has the minimum synchronization cost if C = 
ai{n — 1) for any (F, H, s) with TZ satisfying (A-1) and (A-II). Thus, the 
optimality is a property that depends solely on the network structure and not 
on the dynamical component of the system. 

Theorem 4 The following statements are equivalent: 

(i) A network has the maximum synchronizability. 

(ii) A network has the minimum synchronization cost. 

(iii) The Laplacian eigenvalues of a network satisfy condition (36) . 



PROOF, (iii) =^ (i) and (iii) (ii) are precisely what Corollary 3 states, so 
we are left with proving (i) =>■ (iii) and (ii) =^ (iii). We do this by showing 
their contrapositivcs, i.e., that if (iii) docs not hold, then neither (i) nor (ii) 
hold. If (iii) does not hold, we either have Re A2 < Re A„ or Im 7^ for 
some k, where 2 < A; < n. If Re A2 < ReA„, then, by (30) and (31), we have 



02 Re A2 CX2 

— ' \ ^ * 

«! Re An «! 



(40) 



and 



C > 



ReA- 



> ai 



(n - 2) + 



Re A„ 
ReA. 



(41) 



> ai{n — 1), 



and hence neither (i) nor (ii) hold. If Re A2 = • • • = Re A„ and Im A^ 7^ for 
some k, then there exist systems with TZ satisfying (A-I) and (A-II) (in fact, 
most systems; see Fig. 3) such that {Re (o-Afc) | aXk G TZ} is an interval with 
endpoints at a[ and where cti < a[ < a'2 < 0.2, implying 



an a2 
-T < — ■ 



(42) 



Here we used the fact that T^reai is an interval with endpoints at and a2, 
as in the proof of Corollary 3. Again, by using the same argument as in the 
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Fig. 3. A situation leading to suboptimal S and C in the proof of Theorem 4. 
proof of Theorem 1 with ai and a2 replaced with a'l and a'2, we have 



so neither (i) nor (ii) hold. We have shown (i) <^ (iii) and (ii) 4^ (iii), which 
prove the equivalence of the three statements. □ 

The conclusions about the synchronization cost C in Theorem 1, Theorem 2, 
Corollary 3, and Theorem 4 remain valid if the assumptions on the stability 
region are weakened, with straightforward modification to the proofs. The 
bound on C in Theorem 1 and the identity for C in Theorem 2 are valid 
without assuming (A-Il). The conclusion about C in Corollary 3 remains valid 
if (A-I) and (A-11) are replaced with the condition that cti (as a point in the 
complex plane) is contained in T^-reai- The latter condition is also sufficient for 
the equivalence between statements (ii) and (iii) in Theorem 4. 

In view of the equivalence in Theorem 4 and for convenience, we will use the 
following terminology in the rest of the paper. 

Definition 5 We say that a network given by L is optimal if it satisfies 
condition (36). 

Thus, the class of optimal networks is the set of simultaneous solutions to 
the two optimization problems. The uniform global coupling, in which all 
oscillators are connected to all the other oscillators with weight X/n on all 
links, is perhaps the simplest example of an optimal network [Fig. 4(a)]. The 
non-zero Laplacian eigenvalues in this case are A2 = • • • = A„ = A, and 




(43) 



and 



C > a'i(n - 1) > ai{n - 1), 



(44) 
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Fig. 4. Simple examples of optimal networks, (a) The global coupling configuration 
with uniform weight A/n on every link, where n = 6 and each double arrow repre- 
sents two links, one in each direction, (b) The outward oriented star configuration 
with uniform weight A on every link, where each arrow represents a single link. 

therefore this network satisfies condition (36). Another simple example is the 
outward oriented star, defined here as having a single node connected to all 
the other nodes with uniform weight A on all of these links [Fig. 4(b)]. The 
Laplacian eigenvalues are also A2 = • • ■ = A„ = A for this configuration. 
However, there are many more networks that are optimal, as we will see in 
Section 6. 



5 Diagonalizability of Optimal Structures 

We now study the diagonalizability of the optimal networks defined by (36). 
Here we say that a network is diagonalizable if the corresponding Laplacian 
matrix is diagonalizable. Otherwise, the network is called non-diagonalizable. 
We will show that all networks that are diagonalizable must satisfy a special 
structural condition. This has a rather surprising consequence that most op- 
timal networks are non-diagonalizable. This means that the extension of the 
master stability analysis in Section 2 was indeed necessary for a proper treat- 
ment of the optimization problems. The following theorem gives a complete 
characterization of the optimal networks that are diagonalizable. 

Theorem 6 The following two statements about a given network are equiva- 
lent: 

(i) The network is optimal and diagonalizable. 

(ii) The oscillators can be divided into two groups: those with uniform positive 
output to all the other oscillators and those with no output at all. In 
addition, there is at least one oscillator in the first group. 

PROOF. Suppose that a network satisfies condition (i), i.e., L is optimal, 
so that A = A2 = • • • = A„ > 0, and L is diagonalizable. Then, the eigenspace 
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associated with the eigenvalue A has the maximum possible dimension of n — 1, 
and so docs the solution space of the eigenvalue equation Lx = Xx, which can 
also be written as (L — XI)x = 0. This implies that all rows of the matrix 
L — XI must be constant multiples of the first row, so the matrix must be of 
the form 

/ „ „ \ 



L- A/- 



C2ai C2a2 



C20'n 



(45) 



\Cnai c„a2 ■ ■ ■ c„a„ J 
Then, the Laplacian matrix itself takes the form 



Oi + A 02 • • • (In 

C2ai C2a2 + A ■ • • C2an 



y c„ai c„a2 • • • c„a„ + Ay 



However, from the property of a Laplacian matrix that J2j Lij 
that C2 — • • ■ — Cn — ^1 and hence 



(46) 



0, it follows 



ai + A a2 
ai a2 + A 



On 

an 



\ 



02 



a„ + A 



(47) 



The off-diagonal entries in the jth column represent the strength of the output 
from the jth oscillator to the other oscillators. Hence, this form of L implies 
that each oscillator j either has the same positive output to all the oscillators 
{ttj 7^ 0) or it has no output at all {aj — 0). In addition, there must be at 
least one oscillator that has positive output, since otherwise the network is 
completely disconnected, and it is impossible for it to be optimal. 

Now suppose that the network satisfies condition (ii). Then, since the strength 
Aij of the connection from node j to node i ^ j depends only on j, the 
adjacency matrix A must have the form 



62 
bi 



bn 



ybi 62 



(48) 
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(b) 




Fig. 5. Examples of optimal networks that are diagonalizable. Each arrow represents 
a single directed link. Thick, medium, and thin arrows have weight A/2, A/3, and 
A/6, respectively. For each network, all nonzero eigenvalues are A. 

where 6^ > 0. Using J2j Lij = 0, we can show that the Laplacian matrix L 
must be in the form 



A - 6i -62 
-bi A - 62 



-br. 
-br, 



\ 



-bo ■ ■ ■ \ — br 



(49) 



where X = J2ibi > 0- From this, it is straightforward to show that the eigen- 
values of L are = Ai<A2 = -- - = An = A, and the eigenspace of A has 
dimension n — 1. Thus, the network is optimal and diagonalizable. □ 



The uniform global coupling and the outward oriented star configurations, 
shown in Fig. 4, are examples of optimal networks that are diagonalizable. 
Figure 5 shows two more such examples. For both examples, the Laplacian 
eigenvalues A2, . . . , Ay are all equal to A and the corresponding eigenspace has 
the full dimension of 6. 

In fact, condition (ii) in Theorem 6 provides us with an explicit procedure 
to construct all optimal networks that are diagonalizable. Start with a set of 
n nodes with no links and choose nonnegative numbers &i, 62, ■ ■ ■ , bn, not all 
zero. For each i, make directed links from node i to the n — 1 other nodes with 
weight bi on each link. Regardless of the choice of 6j's, the resulting network is 
guaranteed to satisfy condition (ii), and different choices of 6j's are guaranteed 
to generate all such networks. In the example of Fig. 5(a), the three nodes in 
the middle have links to all the other nodes with weight bi — A/3, while the 
remaining nodes have no outgoing links {bi = 0). In Fig. 5(b), the three middle 
nodes have outgoing links with different weights {bi = A/2, A/3, A/6) that add 
up to A. In this class of networks, the {n — l)-degenerate eigenvalue is always 
equal to J2i bi- 
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The theorem imphes that for any optimal network that is diagonahzable, there 
is at least one oscillator with uniform positive output to all the other oscil- 
lators in the network. This, however, is unlikely to occur in a large realistic 
complex networks, implying that if synchronization is important for such net- 
works, they are hkely to be non-diagonahzable, or at least close to being 
non-diagonalizable. 



6 Optimality of Hierctrchical Network Structures 

Here we present another subclass of optimal networks, characterized by three 
structural conditions that are more intuitive than (36). The first condition 
ensures connectedness of the network, the second ensures well-defined hierar- 
chy of nodes, and the third ensures uniformity of total input strength in each 
node. To conveniently state the first condition, we define an oriented spanning 
tree to be a directed subnetwork that is a tree and spans the entire set of 
nodes, with the links oriented in the direction from the root to the branches 
of the tree. Thus, the existence of an oriented spanning tree embedded in a 
network is equivalent to the existence of a node from which all other nodes 
can be reached by following the directed finks. 

Theorem 7 Suppose that a network satisfies the following three conditions: 

(i) It embeds an oriented spanning tree. 

(ii) It has no feedback loops. 

(iii) For all nodes that receive positive input, the sum of input strength YliO^i — 
Ljj is equal to a constant X. 

Then, the network is optimal and the {n — l)-degenerate eigenvalue is equal to 
X. 



PROOF. Suppose that a network satisfies the given conditions (i)-(iii). By 
starting from an arbitrary node and traversing nodes following links in the 
reverse direction, we must eventually either return to a node already visited, 
thus creating a feedback loop, or arrive at a node without any input. By 
condition (ii) wc cannot have any feedback loops, so we must arrive at a node 
that receives no input. Such a node can only be the unique node at the root of 
the oriented spanning tree that is guaranteed to exist by condition (i), since 
any other nodes in the tree must be reachable from the root. 

Now let us assign the index 1 to the unique node without input. Consider the 
network obtained by removing the node 1 and any links from it. Applying the 
same argument as above to this subnetwork (but now only the existence part) , 
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Fig. 6. (Color online) (a) Example of optimal network with hierarchical structure. 
Thick, medium, and thin arrows have weight A, 2A/3, and A/3, respectively, where 
the sum of input strengths in each node is normalized to A. The nodes are numbered 
according to the ranking and colored by the levels in the hierarchy, (b) The nodes 
in (a) are rearranged to make the hierarchical levels more clearly visible. 

we see that there is at least one node without input within the subnetwork. 
The only input to such a node is from node 1. Let n2 be the number of 
such nodes. We arbitrarily index them 2, 3, . . . , 71,2 + 1. Let us then consider 
the network obtained by removing these n2 additional nodes, together with 
all associated links. Applying the same argument again, we see that there 
is at least one node whose only input is from nodes 1, 2, . . . , 77,2 + 1, and we 
index them n2 + 2,n2 + 3, . . . ,n2 + + 1. Repeating this argument until we 
assign indices to all the nodes, we obtain an indexing in which all the links 
are from a node with smaller index to a node with larger index. This means 
that the Laplacian matrix L of the network using these indices is a lower 
triangular matrix, and hence the diagonal elements are its eigenvalues. Since 
the diagonal elements of L are precisely the total input strength of the nodes, 
the first one is 0, corresponding to the unique node without input, and all the 
others are A, which follows from condition (iii). Hence the network satisfies 
the condition (36) and therefore is optimal, by Theorem 4. □ 

Note that condition (i) of Theorem 7 is equivalent to Re A2 > 0, which follows 
immediately from a recent result in Ref. [43], and this generalizes the notion 
of connectedness to directed networks. Condition (i) is necessary for a network 
to satisfy the stability condition (10). In other words, the network must be 
connected in this sense to make sure that it is at least compatible with the 
possibility of stable complete synchronization. 

In Fig. 6, we show an example of a network satisfying conditions (i)-(iii) of 
Theorem 7. The method of assigning indices to the nodes described in the 
proof above defines a ranking of the nodes in such a network. In this ranking, 
links exist only from a node of higher rank to a node of lower rank. In addition, 
the method defines a hierarchical structure with multiple levels. The top level 
contains the unique node without input. The second level consists of 77.2 nodes 
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that receive input only from the top IcveL The third level consists of ns nodes 
that receive input only from the top two levels. The rest of the hierarchical 
structure is defined similarly, so that links exist only from a higher level to a 
lower level. These hierarchical levels are indicated by different colors in Fig. 6, 
and in panel (b), the nodes are rearranged to make the levels more clearly 
visible. 

The flow of information about the dynamical state of the oscillators in this 
hierarchical structure is unidirectional; it flows only from the top down to the 
bottom. With this picture in mind, the reason for the optimality of such a net- 
work can be understood intuitively as follows. The top node in the hierarchy 
receives no input and acts as the "master" oscillator that dominates the net- 
work dynamics. If the coupling strength a is chosen so that A((tA) < 0, where 
A > is the (n — l)-degenerate eigenvalue, then the oscillators in the second 
level, which receive input only from the master, will synchronize themselves 
with the master. An oscillator in the third level, which receive input only from 
those in the top two levels, must also synchronize, since normalization of the 
total input strength makes the equation effectively look as if it were receiv- 
ing input from a single oscillator synchronized with the master. Repeating the 
same argument for the rest of the hierarchical levels in the network, we see that 
all oscillators must eventually synchronize. Thus, conditions (i)-(iii) guaran- 
tee stable synchronization in the entire range of a such that A((tA) < 0. This 
makes perfect sense because the stabihty condition (10) becomes A((tA) < 
when the optimality condition (36) is satisfied. 

Notice that this argument is very similar to the argument in Section 2 that 
was used to derive the stability condition for non-diagonalizable cases. This 
suggests that the networks satisfying conditions (i)-(iii) may also suffer from 
long transient before converging to the synchronized state. For these networks, 
the number of levels in the hierarchy is strongly related to the length of the 
transient. In addition, the similarity suggests that most of these hierarchical 
networks in Theorem 7 are non-diagonalizable. In fact, the only exception is 
when there are only two levels in the hierarchy, with one top node connected 
to all the other nodes with uniform weights, leading to the outward oriented 
star configuration. 

Theorem 8 Let L represent a network satisfying (i)-(iii) in Theorem 7. Then, 
L is diagonalizable if and only if it is the outward oriented star. 



PROOF. Suppose L is diagonahzable. Then, by Theorem 6, for each i, oscil- 
lator i either has the same nonzero output to all the other oscillators or has no 
output at all. From the argument in the proof of Theorem 7, there is a unique 
oscillator without any input. This oscillator must have uniform output to all 
the other oscillators, since it would be isolated otherwise. Any of the other 
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oscillators must have no output at all, since otherwise it would have output to 
all the other oscillators, including the first one, which leads to a feedback loop. 
Thus, we have have the outward oriented star configuration, and the weights 
on the links are uniform because of the condition (iii) in Theorem 7. 

If the network is the outward oriented star, then it is clear that each oscillator 
either has the same nonzero output to all the other oscillators or has no output 
at all. Therefore, by Theorem 6, it must be diagonalizable. □ 



This result is intimately related to the structure of branches in the underlying 

spanning trees. In an oriented spanning tree, it can be shown that the number 
of independent eigenvectors associated with an (n — l)-degenerate eigenvalue 
A > is equal to the number of branches in the tree. 



7 Optimization with Topological Constraints 

Let us now consider optimization problems with topological constraints. Sup- 
pose that the oscillators arc constrained to interact only within a given network 
topology represented by a symmetric matrix Aq defined by 



Note that Aq represents an undirected network of interaction topology. To 
make the system compatible with synchronization, we assume that this net- 
work is connected, i.e., there is a path between any two nodes. The problem is 
to choose an assignment of weights and directions for the hnks in this network, 
so that the resulting network is optimal. Let Wij > be the weight assigned 
to the directed link from j to i. With this assignment, we obtain a network 
with adjacency matrix A given by Aij = (y4o)ijW^ii- Then, the constrained 
optimization problems can be formulated as follows. Given a connected net- 
work topology ^0 of allowed interactions, for which assignment W does the 
resulting network have the maximum synchronizability and/or the minimum 
synchronization cost? 

Remarkably, there always exists an assignment that achieves the optimality 
defined by (36) for any constraining topology Aq that is connected. We can 
explicitly construct solutions using Theorem 7, together with the fact that we 
can always find an oriented spanning tree within the topology Aq. Indeed, if 
we properly assign directions to links along such a tree, then conditions (i) 
and (ii) are clearly satisfied. The fact that properly weighted trees can have 




1, if distinct oscillators i and j are allowed to interact, 
0, otherwise. 



(50) 
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Fig. 7. (Color online) (a) Example of weight and direction assignment within a 
given interaction topology based on an oriented spanning tree constructed by the 
breadth- first search. Arrows indicate directed links with nonzero weight, while 
dashed lines are links with zero weight, (b) The nodes in (a) are rearranged to 
make the hierarchical levels more clearly visible. 

identical nonzero eigenvalues has been noted before, without considering their 
non-diagonahzability [44]. 

One way to explicitly construct an oriented spanning tree is the well-known 
procedure called the breadth-first search. The procedure also determines a 
ranking of the nodes and the hierarchical levels. First, we choose an arbitrary 
node as the top rank node which forms the top hierarchical level by itself. 
Then, we find all nodes that can be reached from the first node, make con- 
nections to them, and rank them arbitrarily following the first node. These 
nodes form the second level in the hierarchy. The third level consists of the 
nodes that are not yet discovered but can be reached by following two links, 
and again we rank them arbitrarily following the nodes already discovered. 
We make connections to these nodes, making sure to choose only one con- 
nection to each node. We repeat this until we discover all the nodes in the 
network, and the resulting directed network is an oriented spanning tree. This 
procedure can produce at least n distinct oriented spanning trees, one for each 
choice of the root node, but in general there are many others, and many of 
them cannot be found by this procedure. Indeed, from the Matrix- Tree The- 
orem it follows that the number of all such oriented trees is H^^g/^j) where 
/i2 , • • • , A'-n are the nonzero Laplacian eigenvalues of the underlying undirected 
network defined by matrix Aq. For a globally connected network, for example, 
the number of oriented spanning tree is n"~^, which is a huge number even 
for relatively small networks. 

To ensure that condition (iii) is also satisfied, we assign the same weight to 
all the links in the oriented spanning tree. The resulting network is guaran- 
teed by Theorem 7 to have the maximum synchronizability and the minimum 
synchronization cost (see Fig. 7 for an example). Therefore, for a given con- 
nected topology of possible interactions, there are at least as many solutions 
to the constrained optimization problems as the number of oriented spanning 
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trees. However, there are certainly many more ways to assign weights and 
directions so as to satisfy the conditions in Theorem 7. For example, any ad- 
dition of directed links to an oriented spanning tree, that is allowed by the 
topology of possible interactions and that does not create loops, leads to an 
optimal network after normalization of inputs. It is rather remarkable that 
these solutions allow the network with an arbitrary topological constraint to 
achieve the same synchronization efficiency — from both synchronizability and 
cost viewpoints — as the global coupling configuration. 

An interesting property of this construction is that the choice of the "master" 
oscillator does not matter in achieving optimality. Despite the intuition that 
the nodes with the largest number of links in the given topology are the most 
natural choices for the master, the above construction shows that any node 
can be the master. Moreover, the direction of each link in an optimal network 
is not necessarily related to the properties of the nodes it connects. 

Because the weight assignments based on oriented spanning trees explore only 
some of the many potential links, the optimality of the resulting networks can 
be interpreted as a synchronization version of the paradox of Braess for traffic 
flow [45,46], in which removing links leads counter- intuitively to improved 
performance of the network. It is interesting to notice that similar directed 
networks without feedback loops also emerge as gradient networks [47] and in 
the study of transport processes on complex networks [48] . 

An immediate consequence of Theorem 6 is that all solutions of the constrained 
optimization problem are non-diagonalizable, unless some oscillator is allowed 
to interact with all the other oscillators. 

Corollary 9 Suppose that the topology Aq of interactions allows no oscilla- 
tor to interact with all the other oscillators. If an assignment of weights and 
directions leads to an optimal network (which is always possible as long as Aq 
is connected), then the resulting network is non-diagonalizable. 

PROOF. Suppose that the optimal network obtained by the assignment of 
weights and directions is diagonahzable. Theorem 6 implies that there must 
be at least one oscillator having nonzero output to all the other oscillators. 
Therefore, under the assumption of the Corollary, the network must be non- 
diagonalizable. □ 

Once again, this shows that the extension of the master stability framework 

was indeed necessary to treat the problem. The global coupling topology and 
star topology are among the exceptional cases where an oscillator can com- 
municate with all the other oscillators, but such a situation is uncommon in 
a large complex network. 
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The oriented spanning trees form a subclass of all optimal networks con- 
strained under the given topology of interactions. One property that sets the 
oriented spanning trees apart from others in the class of optimal networks is 
that in addition to maximizing synchronizability and minimizing the synchro- 
nization cost, it also minimizes the number of links used {i.e., the number of 
links with nonzero weight). This is because the number of links in an oriented 
spanning tree, which is always n — 1, is the minimum number of links required 
to span all n nodes. 

The oriented spanning trees can be used as a basis for optimization under even 
stricter constraints. Suppose that, in addition to constraining the topology of 
interactions, we have constraints on the directions of allowed links. In other 
words, we consider optimization among all subnetworks of a given directed 
network represented by 



Note that in this case, Aq is not necessarily symmetric. As long as the directed 
network given by is connected in the sense that it embeds an oriented span- 
ning tree, explicit construction of an optimal subnetwork is possible. These 
optimal subnetworks include the embedded oriented spanning tree itself and 
any other subnetwork satisfying the conditions in Theorem 7. 

Interestingly enough, despite all the optimal properties that stem from the 
properties of oriented spanning trees, undirected tree topology was found to 
be among the most difficult to synchronize [49] . Moreover, the out-degree dis- 
tribution of these optimal networks can be highly heterogeneous, in sharp con- 
trast with the case of undirected networks [24] . These highhght the significance 
of directionality of the interactions in determining the synchronizability of net- 
works. The fact that directed networks may have advantage over undirected 
ones is consistent with the finding in [27-29,50] that asymmetric coupling in 
networks of chaotic oscillators has positive effect on synchronization. 



8 Concluding Remeirks 

In this work, we have considered the problem of maximizing the synchroniz- 
ability and minimizing the synchronization cost of oscillator networks. By ex- 
tending the master stability formalism to the case of non-diagonalizable Lapla- 
cian matrices, we have shown that the solutions of the optimization problems 
can be completely characterized by the simple condition (36) on the Lapla- 
cian eigenvalues. We have also shown that the intuitive structural conditions 
(i)-(iii) in Theorem 7, which facilitate unidirectional information flow and nor- 




1, if oscillator i may receive connection from oscillator j i, 
0, otherwise. 



(51) 
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(a) (b) 





(d) 



Fig. 8. Examples of optimal networks that are neither diagonalizable nor hierarchi- 
cal. An arrow indicates a directed link, while a double arrow indicates two directed 
links. Thick, medium, and thin arrows have weight A, A/2, and A/3, respectively. 
For each network, all nonzero eigenvalues are A. 

malized input strength, guarantee optimality. In addition, we have considered 
optimization under topological constraints and shown that we can explicitly 
construct optimal networks using oriented spanning trees. Furthermore, by a 
complete characterization of diagonalizable optimal networks, we have proved 
that most optimal networks are actually non-diagonalizable, which necessi- 
tates the extension of the master stability formalism. Since spectral analyses 
are also relevant for other dynamical processes on networks (e.g., diffusion and 
other spreading processes), it would be interesting to see these results applied 
to the study of different network phenomena. 

Structural properties of optimal networks, such as those given in our theo- 
rems for the diagonalizable and hierarchical networks, can serve as a general 
guideline for designing networks for synchronization. For such applications, 
however, it is important to address the question of robustness. What is the 
effect of structural perturbations on the optimality of these networks? We ex- 
pect that a perturbation theory can be used to show that small deviations from 
the optimal structures will induce only a small change in S and C. Another 
type of robustness question is about the effect of dynamical perturbations: 
how does disturbances introduced at one or more nodes in the synchronized 
state propagate over the network? 

Within the class of all optimal networks, we have explicitly constructed a large 
number of important networks: those that are diagonalizable (Theorem 6) and 
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Fig. 9. Example of a suboptimal network resulting from combining two optimal 
networks. The lower-level node of a hierarchical network with two nodes is replaced 
by a global coupling network with two nodes. We can show that no combinations 
of weights on the two downward links can make this network optimal. 

those that are hierarchical (Theorem 7). However, there are many other opti- 
mal networks that fall into neither categories (sec Fig. 8 for examples). Some 
of them can be constructed by combining two optimal networks. For example, 
the network shown in Fig. 8(d) is a combination of a hierarchical network 
and a diagonalizable network, constructed by replacing the root node in the 
5-node outward oriented star (hierarchical) with the 3-node global coupling 
configuration (diagonalizable). However, the other examples in Fig. 8 cannot 
be constructed in a similar fashion. There are also cases in which this kind 
of construction does not lead to an optimal network (see Fig. 9 for an exam- 
ple). The explicit construction of all optimal networks is an important open 
problem. 

The complete characterization of the entire class of optimal networks is ex- 
pected to have a profound implication on the widely assumed hypothesis that 
synchronizability plays an important role in the evolution of many real-world 
complex networks. If the signatures typical of optimal networks are found 
to be absent in these networks, then the hypothesis may be questionable; 
synchronization may be less significant than other competing factors such as 
robustness, or else, the oscillator network model does not describe the essen- 
tial features of the system. Hierarchical structures have been suggested to play 
a significant role in the network of motor neurons of Aplysia [51] and in the 
mechanism of invariant visual representation in the human brain [52]. Hier- 
archical structures are also common in many human organizations, perhaps 
because they can better facilitate coordinated activities. It is important to 
examine the existing real data as a first step toward testing this and possibly 
other hypotheses about the evolution of complex networks. 
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